


#load packages:
library(R.utils)
library(raster)
library(sf)
library(e1071)

#set project name:
projectName<-"WHOSnakes"

#define directories:
MyDir<-paste("/home/jc217070/Maxent",sep="")
LayersDir<-paste("/home/jc217070/Maxent/Layers")
SDMDir<-paste(MyDir,"/Outputs3_Final/Final_Individual_Spp",sep="")
SWDDir<-paste(MyDir, "/Maxent_SWDs/spp_occs",sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]

worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries_2024.shp",sep=""))  #maptools
mycntrs<-c("SOUTH SUDAN")
#CountryPolies<- subset(worldmap, worldmap$NAME_CTY %in% mycntrs)     #get subset of countries that species occurs in
CountryPolies<- subset(worldmap, worldmap$ADM0_NAME %in% mycntrs)     #get subset of countries that species occurs in
myext<-extent(CountryPolies)

#extract records for all species for this region:
CountryRecs<-data.frame()
for (sp in spp){
  
  print(paste(Sys.time(),"extracting records for",sp,":"))
  
  n<-which(spp==sp) 
  print(paste(Sys.time(),"extracting country records for", n, sp, sep=" "))
  
  if(file.exists(paste(SWDDir, "/",sp,"_allRecs.csv",sep=""))){
  sprecs<-read.table(file=paste(SWDDir, "/",sp,"_allRecs.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
  coords<-data.frame(long=sprecs$long,lat=sprecs$lat)  # get coordinates of occurrence records
  spcoords <- SpatialPoints(coords[,1:2],proj4string=CRS("+init=epsg:4326"))
  
  spcoords2<-st_as_sf(spcoords)
  
  #clip records to country
  addrecs<-st_filter(spcoords2, CountryPolies)
  addrecs<-as.data.frame(st_coordinates(addrecs))
  names(addrecs)<-c("long","lat")
  
  if(length(addrecs[,1])>0){
    addrecs$unloc<-paste(round(addrecs$lat,1),";",round(addrecs$long,1),sep="")
    finrecs<-sprecs[sprecs$unloc %in% addrecs$unloc,]
    CountryRecs<-rbind(CountryRecs,finrecs)
  }
  
  write.table(x=CountryRecs, file=paste(MyDir,"/SouthSudanClip/", gsub(" ","-",paste(mycntrs,collapse="_")),"_allRecs.csv",sep=""),
              qmethod = "double", na = "NA", row.names = FALSE, col.names = TRUE, sep=",",quote=TRUE)
  
  }
}

#get my region polygons and rasterize
GLOBraster<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")
snapraster<-crop(GLOBraster, myext)
values(snapraster)<-ifelse(values(snapraster)!=-Inf,0,-Inf) #raster showing all land as 0, ocean as NA
snapraster[snapraster == 0] <- NA
regionras  <- rasterize(x=CountryPolies,y=snapraster,as.numeric(1))

#find species with polygons overlapping this region (according to experts, not SDMs)
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
sf_use_s2(FALSE)
mypoliesTF<-st_intersects(WHOpolies,CountryPolies,sparse=FALSE) #species polygons overlapping with this region polygon
WHOpolies$TF<-mypoliesTF
mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region

#extract SDMs for all species for this region:
for (sp in gsub(" ","_",myspp[16:length(myspp)])){
  
  print(paste(Sys.time(),"extracting rasters for",sp,":"))
  
  outfile<-paste(SDMDir, "/",sp,"_Current_CropThreshCDCut.asc",sep="")
  outfile1<-paste(SDMDir, "/",sp,"_Current_Risk.asc",sep="")
  outfile2<-paste(SDMDir, "/",sp,"_Current_Exp.asc",sep="")
  outfile3<-paste(SDMDir, "/",sp,"_Current_Imp.asc",sep="")
  
  if(file.exists(paste(outfile,".gz",sep=""))){
    
    gunzip(filename=paste(outfile,".gz", sep=""),
           destname=paste(outfile),
           remove=FALSE, overwrite=TRUE)
    
    SDM<-raster(paste(outfile), crs="+init=epsg:4326")
    SDM<-SDM*regionras
    writeRaster(SDM, filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_SDM.asc",sep=""), format="ascii", overwrite=TRUE)
    gzip(filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_SDM.asc",sep=""))
    
   ###
    
    gunzip(filename=paste(outfile1,".gz", sep=""),
           destname=paste(outfile1),
           remove=FALSE, overwrite=TRUE)
    
    Risk<-raster(paste(outfile1), crs="+init=epsg:4326")
    Risk<-Risk*regionras
    writeRaster(Risk, filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Risk.asc",sep=""), format="ascii", overwrite=TRUE)
    gzip(filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Risk.asc",sep=""))
    
    ###
    
    gunzip(filename=paste(outfile2,".gz", sep=""),
           destname=paste(outfile2),
           remove=FALSE, overwrite=TRUE)
    
    Exp<-raster(paste(outfile2), crs="+init=epsg:4326")
    Exp<-Exp*regionras
    writeRaster(Exp, filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Exp.asc",sep=""), format="ascii", overwrite=TRUE)
    gzip(filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Exp.asc",sep=""))
    
    ###
    
    gunzip(filename=paste(outfile3,".gz", sep=""),
           destname=paste(outfile3),
           remove=FALSE, overwrite=TRUE)
    
    Imp<-raster(paste(outfile3), crs="+init=epsg:4326")
    Imp<-Imp*regionras
    writeRaster(Imp, filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Imp.asc",sep=""), format="ascii", overwrite=TRUE)
    gzip(filename=paste(MyDir,"/SouthSudanClip/",sp,"_Current_Imp.asc",sep=""))
    
    if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
      file.remove(outfile)
    }
    if(file.exists(paste(outfile1,".gz", sep="")) & file.exists(outfile1)){
      file.remove(outfile1)
    }
    if(file.exists(paste(outfile2,".gz", sep="")) & file.exists(outfile2)){
      file.remove(outfile2)
    }
    if(file.exists(paste(outfile3,".gz", sep="")) & file.exists(outfile3)){
      file.remove(outfile3)
    }
    
  }
}


